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Abstract 

in ' 

^■s^ ' In networks of periodically firing spiking neurons that are interconnected with chemical 

synapses, we analyze cluster state, where an ensemble of neurons are subdivided into a few 

1 i ■ clusters, in each of which neurons exhibit perfect synchronization. To clarify stability of clus- 

ter state, we decompose linear stability of the solution into two types of stabilities: stability 
Qh , of mean state and stabilities of clusters. Computing Floquet matrices for these stabilities, we 

r*j ' clarify the total stability of cluster state for any types of neurons and any strength of inter- 

* '""j ' actions even if the size of networks is infinitely large. First, we apply this stability analysis 
C"| . to investigating synchronization in the large ensemble of integrate-and-fire (IF) neurons. In 

one-cluster state we find the change of stability of a cluster, which elucidates that in-phase 
synchronization of IF neurons occurs with only inhibitory synapses. Then, we investigate 
entrainment of two clusters of IF neurons with different excitability. IF neurons with fast 
decaying synapses show the low entrainment capability, which is explained by a pitchfork bi- 
\f~^ , furcation appearing in two-cluster state with change of synapse decay time constant. Second, 

f*"^ ■ we analyze one-cluster state of Hodgkin-Huxley (HH) neurons and discuss the difference in 

' synchronization properties between IF neurons and HH neurons. 

o : 

m ■ 

1 Introduction 

It has been revealed that periodically firing interneurons exhibit in-phase synchronization during 
the gamma oscillations (20-80 Hz) and the sharp wave burst (100-200 Hz) [1]. Interneurons are 
found to be connected through inhibitory chemical synapses. Therefore, a significant effort has 

• *h . been devoted to understand a role of inhibitory chemical synapses in in-phase synchronization in 
/\ ' a large ensemble of neurons [2]. One major analytical approach to investigate synchronization 

of neurons is the phase reduction method, in which behavior of periodically firing neurons are 
reduced to the simple phase dynamics [3-6]. This phase reduction method is, however, applicable 
only to the case of weak couplings. To understand a role of strong couplings in synchronization 
of neurons we have to adopt different approach. 

One difficulty in investigating strongly coupled neurons is time delayed interactions due to 
chemical synapses. Taking account of these time delayed interactions Hansel et al. have com- 
puted Floquet matrix and analyzed synchronization in a couple of strongly coupled neurons [4]. 
The size of this Floquet matrix, however, increases as the size of neural networks increases. 
Therefore, it is difficult to apply their approach to investigating the large size of neural networks. 

Bressloff et al. have presented another scheme to deal with chemical synapses, which al- 
lows us to analyze stability of networks of integrate-and-fire (IF) neurons without computing the 
explicit form of Floquet matrix [5]. In some large size of neural networks they have found the 
degeneracy of eigenvalues, which makes it easy to analyze synchronization of many IF neurons. 
Actually, such degenerate eigenvalues in stability analysis are found not only in IF neurons but 
also in many synchronization phenomena induced by mean field interactions. A most prominent 
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example of this degeneracy is seen in synchronization in an ensemble of chaotic oscillators such 
like Lorenz equations and logistic maps [7-10]. Just using the general properties of mean field 
interactions we can decompose linear stability of synchronous state of chaotic oscillators to two 
different components, which define the so called tangential Lyapunov exponents and transversal 
Lyapunov exponents. It must be noted that the result of this decomposition clearly indicates the 
occurrence of degeneracy regarding transversal Lyapunov exponents. Synchronization in many 
chaotic oscillators is thus characterized by only a small number of exponents included in tangen- 
tial and transversal Lyapunov spectrum even if the system size is infinitely large. 

In the present paper we employ these sophisticated reduction techniques in the chaos syn- 
chronization theory to investigate synchronization in the large number of neurons. The target of 
the analysis is cluster state, where an ensemble of neurons are subdivided into a few clusters, in 
each of which neurons exhibit perfect synchronization. To evaluate the degeneracy of eigenval- 
ues we carry out the above-mentioned decomposition of a linear stability and define stability of 
mean state (tangential Floquet multipliers) and stabilities of clusters (transversal Floquet multi- 
pliers). Stability of mean state elucidates if cluster state is stable in the dynamics among clusters 
while stabilities of clusters clarify whether small perturbations in each cluster converge to van- 
ish. We explicitly compute Floquet matrices of these stabilities for arbitrary neuron dynamics. 
Therefore, we can elucidate stability of cluster state for any types of neurons, even if the size of 
networks is infinitely large and neurons are connected through strong couplings. 

To give a concrete example of the present stability analysis we first analyze networks of IF 
neurons interacting through uniform chemical synapses. In this analysis, we find the change of 
stability of a cluster, which elucidates that in-phase synchronization of a large ensemble of IF 
neurons occurs with only inhibitory chemical synapses. In addition, we investigate two clusters 
of neurons with different excitability, and discuss the relationship of their entrainment properties 
to the synapse decay time constant. Second, we analyze one-cluster state of Hodgkin-Huxley 
(HH) neurons and discuss the difference in synchronization condition between IF neurons and 
HH neurons. 

The paper is organized as follows. In Sec.|2]we present the dynamics of neural networks that 
include Q clusters of spiking neurons. In Sec. [3] we present the analysis for cluster state of the 
neural networks. This analysis is applied to networks of IF neurons in Sec.|4] Then, we analyze 
synchronization of HH neurons in Sec. [5] Finally, in Sec. [6] we give a brief summary and discuss 
the future problems that can be solved by the present approach. 



2 Networks of spiking neurons coupled with chemical synapses 

We consider a spiking neuron whose state is represented by 71-dimensional vector 

x= (v,wi,w 2 ,...,w„_i) T , (1) 

where v represent the membrane potential and {w/} /=1 j describe gating of ion channels. 
Typically, the dynamics of a spiking neuron is defined by Hodgkin-Huxley (HH) equations, 
FitzHugh-Nagumo (FN) equations, and so on. We simply represent these neuron dynamics by 

x = F(x). (2) 

In the analysis in Sec. |3j we assume spiking neurons in the form of Eq. (|2ji. Nevertheless, in 
Sec.|4] we will investigate integrate-and-fire (IF) neurons, which cannot be expressed by Eq. (0 
since v of IF neuron changes discontinuously. This discontinuity of IF neuron requires a minor 
corrections of the analysis in Sec. [3] We will discuss this minor correction in Sec.|4] 

We assume that N spiking neurons {x,} are interconnected through chemical synapses. To 
describe the dynamics of synaptic electric currents, we define spike timing by the time when 
membrane potential V; = [x,] i (the first element of vector x,) exceeds the threshold value = 0. 
We represent A:-th spike timing of neurons i by ?, (£), which satisfies 

vi[t i (k)] = [x i [tt{k)]] 1 =6 (3) 

and 

Vi[t i (k)] = [x i [t i (k)]] 1 >0. (4) 
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Then, the dynamics of networks of spiking neurons is expressed as 

x i = F i (x i ) + (I i ,0,...,0y, (5) 

where function F, (x, ) represents the dynamics of neuron i. Variable /, represents a sum of synap- 
tic electric currents, which is defined by 

j=l k=-oo 

where Jjj represents synaptic coupling from neuron j to neuron i, and function S(t) describes 
time evolution of synaptic electric current. We assume S(t ) taking the form 



t <0, 

S ^ = < -J—(e-^-e-^) 0<t. (7) 

Ti - T 2 V / 

where < T2 < T 1 ■ Constants Ti and T2 are termed decay time and rise time, respectively. 

2.1 Neural networks composed of Q clusters of neurons 

In some problems, we have to consider neural networks including several clusters of neurons, 
such like networks including both interneurons and pyramidal neurons. Moreover, we will later 
study entrainment of two clusters of IF neurons that have different excitability between clusters. 
In the present study we analyze neural networks that are composed Q clusters of neurons. We 
assume that neurons share the same dynamical properties within each cluster, that is, we assume 

F,-(x)=F ? (x), i€U q ,l<q<Q, (8) 

where U q represents the set of neurons that belong to cluster q. In addition, we assume that 
synaptic couplings Jii depend only on cluster indexes of pre and postsynaptic neurons, that is, 
we assume synaptic coupling Jn of the form 

Jij=J qq i/N, ieu q ,jeu q ,. (9) 

Substituting Eqs. ^ and (|9j into Eqs. l|3 and we obtain the dynamics of Q clusters of 
neurons: 

Xi = ¥ q (xi) + (I q ,0,...,oy, (10) 

1 2 

h = 77 £ E V L s i f - i e U q . (11) 

JV q '=i. jeUj k 

Note that synaptic electric current in Eq. M li depends only on cluster index q because of the 
assumption in Eq. (|9j. 

3 Analysis 

3.1 Cluster synchronization of periodically firing neurons 

In the present analysis we focus on cluster state, in which spike timing of neurons are written in 
the form 

tt(k) = t*(k)=t*+kr, 

0<t* q <T,ieU q ,q=l,...,Q, (12) 

where asterisks indicates the quantity in stationary state. In this state, neurons emit periodic 
spikes synchronously within each cluster. We further assume that in cluster state not only spike 
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timing but also neuron states are synchronized within each cluster (i.e., x* = x* (i £ U q )). Sub- 
stituting Eq. (1121 into Eqs. JlOb and ( II Q , we obtain the dynamics of stationary state as 

x; = F 9 (x;) + (/,*,0,...,0) T , (13) 

where S(t) is defined by S(t) = JjtS(f + &T) and r 9 = N q /N represents the ratio of the number 
of neurons in cluster q to the total number of neurons. 

To obtain the explicit form of cluster state we have to calculate T and f*,f|, . . . ,t* so as to 
obtain I* and x* It is obvious that we can safely assume t* = 0, and we can calculate remaining 
Q unknown variables: T,t\,t\,...,t*Q from Eqs. Jl 3i and 1141 following the same scheme as our 
previous study [11, 12]. Note that we can compute these variables not only for IF neurons but 
also for general neuron dynamics, as far as the stable cluster state is concerned. 

3.2 Decomposition of linear stability 

To investigate linear stability of cluster state we assume the infinitesimal deviations of state of 
neurons: 

x,= X ; + 5x ; , ieU c/ (15) 
and infinitesimal deviations of spike timing: 

t,{k) = f q {k) + 8t l {k), ieu q . (16) 

From Eq. Q, we obtain 

8n(k) = -8vi[f q (k)}/c q = -[8x i [t* q (k)]] l /c q , 

i € U q (17) 

with 

c q = v^(k)] = [x*[t* g (k)}] l . (18) 

Note that constant c q is independent of k because of the periodicity of the solution. To obtain the 
relation in Eq. dl7> . we must assume continuous neuron dynamics such as HH neurons and FN 
neurons. Note that we have to carry out the more careful calculation in discontinuous dynamics 
like IF neuron as we will discuss in Sec. [4] Expanding the dynamics in Eqs. dlOi and (lilt to the 
first order we obtain the dynamics for the deviations: 

5x,=f;(x;)5x ! + (5/„0,...,0) t , (19) 



8l i = -4l E J q ^S'\t-t* q ,(k)\8tj(kl (20) 



N q< jeUj k 

where F^(x*) denotes Jacobi matrix. 

The naive evaluation of this N x n-dimensional dynamics yields an eigenvalue problem of 
the large size of matrix. Therefore, for each cluster, we define mean state of neurons: 

JV< ? ieu q 

and mean spike timing: 

F #) = TT L *(*)■ ( 22 ) 

Noting Eqs. J19I . ( I20t . and d!7l >. we can write the dynamics for 8x q and 8t q {k) in the closed 
form 

8i q = ¥' q (x* q )8x q + (8l q ,0,---,oy, (23) 
Sl q = -LVvlM' -£(*)! 8 h'( k ) (24) 



'i' 



k 
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with 

St g (k) = -8v q [t* q (k)]/c q = - [Sx^k^/cq. (25) 

Eqs. d23i - d25l define the decomposed stability of the original A^-body stability. We term this 
decomposed stability stability of mean state. It must be noted that stability of mean state in 
Eqs. d23i - (l25t is effectively a problem in a network of Q neurons with couplings J qq tr q i since, to 
the first order, Eqs. J23b - d25l > are equivalent to 

j t (x* + 8x q ) = F q (x* q + 8x q ) + (I q , 0, . . . , 0) T , (26) 

k = L v v L s [ f - 1 « (*) - 5 V (*)] ■ ( 2? ) 

q> k 

Stability of mean state is a necessary condition for the full stability, but not a sufficient con- 
dition. To investigate synchronization of neurons in each cluster we introduce deviations around 
the averaged state: 

x i = x q + 8x i , ieU q . (28) 
Subtracting Eq. d23l from Eq. d!9l > we obtain the dynamics of 5x, as 

Sx i =F' q (x q )8x i , ieU q . (29) 

Eq d29i defines another decomposed stability. We term this decomposed stability stability of a 
cluster. Stability of cluster q is satisfied when N q deviations 8x.j (i G U q ) converge into 0. These 
Nq dynamics are, however, identical. Therefore, it suffices to evaluate one set of deviations to 
determine the stability of one cluster. Note that the determination of the stability of a cluster 
is effectively a problem of a single neuron dynamics under the unperturbed synaptic electric 
current /* since, to the first order, Eq. d29l > is equivalent to 

^ (x* q + 8%) = F q (x* q + 8%) + (/*, 0, . . . , 0) T , i S U q . (30) 



3.3 Floquet matrices for stabilities of clusters 

We can determine stabilities of clusters following the ordinary procedure of Floquet theory. Since 
the solution x* is periodic, F^(x*) is also periodic. Therefore, a solution of Eq. d29i is written in 
the form 

8x i [t: i (k+l)}=M q L 8xt[t* q (k)i ieU q . (31) 

Calculating Eq. d30l > with small initial perturbations we can obtain every elements in matrix . 

n x n matrix M,j- has n eigenvalues {^y } ■ When cluster q is stable, 5x,- (z 6 U q ) must 

converge to zero after a long time. Therefore, the stability of cluster q is fulfilled when the largest 
absolute eigenvalue satisfies the condition 

|\^| <1. (32) 



3.4 Floquet matrix for stability of mean state 

Determination of the stability of mean state is not an easy problem since the calculation of 8l q in 
Eq. G3J 

requires long past deviations of spike timing. To solve this problem, following Hansel et 
al. [4], we introduce the variables: 



Zqi = E4'V E s f -^'( fe ')- 5 v( fc/ ) 



t *,(k')<t 



Zq2 



(33) 
(34) 



«J(*0<» 
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By means of these variables we can exactly rewrite I q in Eq. J27> in the truncated form 



EVv E s t-t;,{h!)-si^) 

q' t*{k)<t*(k>) 



Therefore, 8l q is written as 



Sla 



+s[t-t;(k)] Zq 2[t;(k)}, t*(k)<t. 



-EVv E s ' 5 W) 

q' t' q (k)<t* rf (k') 

-e-^M^Szq^k)} 
-S[t-f q {k)]8z q2 \t* q {k)}, t* q (k)<t. 



(35) 



(36) 



This means that once we know 8z q \[t*(k)] and 8z q 2[t q (k)], we can neglect past deviations of 
spike timing 8t q i(k') that arose before t = t*(k). 

To take the advantage of z q \ and z q 2, we define the vector 



Yq = ([Xq]l,---,[Xq]n,Zql,Z4iy 



(37) 



We safely assume t* < t* +l (q=\,...,Q — 1). Then, since Eq. (124b is equivalent to Eq. ( I36l l. 
from Eqs. d23i . d36l >. and d25l we can show that deviation 8x q [t*(k + 1)] is determined from only 
Sy q [tq{k)] and 8tg,(k gql ) with 



/i q <q' , 
k+l q'<q. 



(38) 



Moreover, deviations 5z^i [t q (k+ 1)] and 5z 9 2[^(^+ 1)] are given as 

8z ql [t* q (k+l)] 

= -E 7 W'V 5 'k( fc + 1) 8J q'( k qq>) 



q 

„-J7 



/*5 V [$(*)] + S(r)5^[f 4 *(*^ 

rEVv^^ ( " +1) ^ (v)]/ ^V(V) 



(39) 



<7 



+ e -^Sz 92 [/*(£)]. 



(40) 



Therefore, we can also determine 8z q \ [t q (k+ 1)] and 8z q 2[t q (k + 1)] from the above mentioned 
variables: 8y q [t*{k)] and 8t q i(k qq >). We can summarize these relationships in the form 



where 



and 



sy q [*,*(*+ 1)] = E V 5 v K> + [£ (*)] . 



A w ' — 



agy^(*+l)] f 1 t () 
d8t q ,{k qq ,) 

dSy q [t* q (k+l)] 

d[8y q [t*(k)]]x 



d8y q [t* q {k+\)\ 
d[8y q [t*(k)]} n+2 



In this equation, A. qq >8y q >[t*,(k qq i)) represents the contribution from 8t q i{k qq i) 



(41) 
(42) 

(43) 
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In Sec. 13.31 we obtain by calculating Eq 



301 with small perturbations. In the sim- 
explicitly for arbitrary neuron dynamics. We obtain 
d8x q [t*(k+ \)]/d8ij{k q< f) smdd8x q [t*'(k+l)}/d[8y g [t*(k)]}i by calculating Eqs. (|26> and 03 



ilar manner, we can compute A q and B 



with small perturbations. Partial derivatives of z q \ and z q 2 have been given in Eqs. d39l > and J40I . 
Therefore, we can compute every elements in matrices A q and B qq i . For the further details of cal 
culation of A q and B qq > see Ref. [12](, though the definitions of A q , B qq >. 
are slightly different from the present ones.) 
We introduce vector 

Y(k) = (y l [t* l (k)Y...y Q [t* Q (k)Yy. 
Then, we can rewrite the relationship in Eq. Hli in the form 



8Y(k+l) =M II SY(&) 



and so on in Ref. [12] 



(44) 



(45) 



with 



where 



/ E 



M 



V 



M" 



E 



*qq 



.M 



^qq+\ 

E 







HQ 



E J 



(46) 



(47) 



Matrix Mg updates 8y q (k) to 8y q (k + 1), and hence matrix updates all the deviations. 
Q(n + 2) x Q(n + 2) matrix has Q(n + 2) eigenvalues {A,"} , in which a triv- 

I ' J l=\,...,Q{n+2) 

ial eigenvalue one is always included as in the case of ordinary Floquet matrix. The stability of 
mean state is satisfied when all other eigenvalues are less than one in absolute value, that is, the 
largest absolute eigenvalue \k\ | and the second largest absolute eigenvalue \X\ | satisfy 



< 1 = A 1 



(48) 



4 Cluster synchronization in networks of integrate-and-fire 
(IF) neurons 

Let us apply the above analysis to networks of IF neurons that are defined as 

Vi = -Vj + v r + I extq + /,, ieU q . (49) 

When V; exceeds the threshold value 0=0, v; is reset to Vo = — 1. The resting potential v r 
is set to v r = 1, which leads intrinsic firing of neurons. We assume that these IF neurons are 
interconnected with uniform couplings: 

*J=f}- (50) 

As we have mentioned, the discontinuity of IF neurons require a minor correction of the sta- 
bility analysis in Sec. [3] Since derivative v, changes discontinuously at spike timing, we define 
c~ = v*[t*(k) —0] and c+ = v*[t*(k) +0]. To take account of discontinuity of v/, we extend the 
perturbed solution v* + 5v,- before/after spike timing f,(£) = t*(k) + 8ti(k) as illustrated in Fig. [2 
and then define 8vj[t*(k)] and 8vf[t*(k)]. These deviations satisfy the condition, 

8 U {k) = -8vr[f q {k)]/c- = -Sv+[t*(k)]/c+. (51) 

We define two types of mean state variables: 8v q = (l/N) Y.ieu q ^ v 7 an< ^ ^ q = (l/N) Y.ieu q $ v t "> 
and two types of deviations ai'ound the mean state: 8vJ = Sv^ + 8\\ and 8vf = 8vf + Svf. 
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Neuron i behaves continuously in the time interval f,(£) < t < ti(k+ 1), during which we can 
carry out the decomposition of linear stability discussed in Sec. [5] Therefore, noting Eqs. j29L 
J49L and i51\ . we obtain 



5vT [/*(*+ 1)] = e-' 8vt\t* q {k)] = 4^' 5vy[ f ; (*)] 



(52) 



Hence, matrix M„ takes the form 



"JL e -T 



From Eq. d32i . we obtain the condition for stability of cluster q: 



'1 „-T 



< 1. 



(53) 



(54) 



Following the similar scheme, we can derive matrices A qcj i and B^. Substituting A qq i and 
into Eqs. J46t and ( I47> yields M", by which we can determine the stability of mean state. 



4.1 One-cluster state of IF neurons (2=1) 

We begin with investigating one-cluster state Q = 1 . In this state, all neurons take part in shaping 
one-cluster in-phase synchronization. One-cluster solution of Eqs. dl 31 and dl4> is found with 
only g < 1 since too strong synaptic couplings with g > 1 leads bursting of neurons. To elucidate 
the stability of the solution with g < 1 , assuming T\ — 3.5, T2 — O.lTj, and /ext.i — 0> we calcu- 
late |Aj|, and |AA| as a function of parameter g as shown in Fig. |2ja). Since the second 

largest absolute eigenvalue of Mil (i.e., |aJ |) is always less than one, the stability of mean state 
is always satisfied. However, the largest absolute eigenvalue of Mf- (i.e., |Aj^|) becomes grater 
than one with excitatory coupling g > 0. Therefore, the stability of a cluster is satisfied with only 
inhibitory coupling g < 0. These results imply that while a self-coupled single neuron (N = 1) 
can exhibit stable periodic firing with both inhibitory and excitatory couplings, in-phase synchro- 
nization of multiple neurons (N > 1) can take place with only inhibitory couplings g < 0. Since 
the networks show the same synchronization properties in all the decay time X\ > 0, Ti — g phase 
diagram takes the simple form as described in Fig.|3b). It turns out that in-phase synchronization 
of a large number of IF neurons occurs with only inhibitory synapses in all the value Ti > 0. 

Figure fJl shows the result of the numerical simulations. While the networks with inhibitory 
couplings (g = —0.5) exhibits the perfect in-phase synchronization, the network with excitatory 
couplings (g — 0.5) settles into the asynchronous state, in which neurons fire periodically with 
uniformly distributed phase shifts. Our stability analysis explains these numerically results well. 



4.2 Two-cluster state of IF neurons (2 = 2 and 4xt,i = 4xt,2 = 0) 

We then investigate two-cluster state Q = 2 for inhibitory coupling g < assuming r\ = = 0.5 
and 7 ext 1 = / ex t 2 = 0. It has been shown that a couple of IF neurons exhibit a pitchfork bifurcation 
with change of synapse decay time constant [4]. We now show that this pitchfork bifurcation 
occurs even in systems of two clusters of neurons. Figure |4] shows Ti — q>2 bifurcation diagram, 
where (fe denotes ti]T . There are three types of solutions: in-phase ((fe = 0,1), anti-phase 
(92 = 0.5), and out-of-phase solutions. Evaluating eigenvalues of M", M^, and M^, we find 
that the solutions denoted by thick lines satisfy the stability of mean state and the stabilities of 
clusters. 



4.3 Entrainment of two clusters of IF neurons with different excitability 

((2 = 2and/ ext)1 =0^/ ext)2 ) 

We extend the above result to investigate the case when the excitability of neurons are different 
between two clusters. Fixing 7 ext ,i = 0, we investigate the behavior of <f>2 with change of 4xt.2 
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in Fig. 13 With 4xt,2 = 0, we find three stable and two unstable solutions, which are consistent 
with the preceding results in Fig.|4] The in-phase solution q>2 = 0, 1 is extended by the change 
of 4xt,2 within the interval — 0.019 </ e xt,2 <0.020. In this interval two clusters of neurons show 
synchronized firing with small phase difference, that is, entrainment occurs. To examine the 
robustness of this entrainment, we plot this range of / ex t,2 as a function of Z\ in Fig. |6] The 
remarkable feature of this phase diagram is the narrow range of 4xt,2 with short decay time 
constant Z\, and it is interesting that the pitchfork bifurcation described in Fig. [4] explains this 
narrow range of / ex t,2- In this bifurcation diagram, the out-of -phase solutions merge into the in- 
phase solutions at Ti = 0. Therefore, the entrained solution in Fig.|5]vanishes in the limit Ti — » 0, 
and this vanishment explains the zero range of 7 ex t.2 at Ti = in Fig. [6] 

On the other hand, the out-of-phase solution (<p2 = 0.5) is considerably robust against the 
change of / e xt2> especially with short T\ (—0.083 </ext,2<0. 080 with i\ = 1.5.) Nevertheless, 
when we apply the external currents to halves of neurons of both clusters (i.e., Q = A,r\ =T2 = 
r 3 = r 4 = 0.25, 4xu = 4xt,3 = 0,4xt,2 = 4xt,4 ¥= °> <Pi = °: <P2 ~ 0, 93 ~ 0.5, <p 4 ~ 0.5), the range 
for successful entrainment is found to be narrow (—0.016 <4xt.2 = 4xt,4 <0.017 with Ti = 1 .5). It 
seems that cluster synchronization easily breaks when we apply heterogeneous external electric 
currents that cause splitting of clusters. 

5 One-cluster state (Q — 1) in networks of Hodgkin-Huxley 
(HH) neurons 

To explore the biological plausibility of synchronization in IF neurons we study more realistic 
neuron model that is defined by HH equations. A HH neuron, whose dynamics is described 
in appendix [X] does not show intrinsic firing without external stimulus. Therefore, we apply 
constant external electric current I ex i = 10 (fiA /cm 2 ) to all of HH neurons and analyze synchro- 
nization in intrinsically firing homogeneous HH neurons assuming the same synaptic couplings 
as Eq. (I50> . Figure shows Ti — g phase diagram, where the condition for stable one-cluster 
state (2=1 and 2 < AO is described. In the large area of inhibitory couplings (g < 0) we find 
stable in-phase synchronization. Beyond Ti =7.0 the behavior of \X\ | and |Afj| is similar to 
those of IF neurons described in Fig.|2ja), and the change of stability occurs at g = because 
of \Xy]\. Below Ti = 7.0, however, synchronization with inhibitory couplings takes place only 
below a certain negative value of g, and excitatory couplings can induce synchronization in some 
conditions. Ti — g phase diagram of IF neurons (Fig. |2jb)) can explain synchronization in HH 
neurons with slowly decaying synapses, though the synchronization condition of HH neurons 
with fast decaying synapses is more complicated than IF neurons. 

6 Discussion 

We have studied cluster state of networks of spiking neurons. We have shown the analytical 
method that can deal with synchronization in the large size of neural networks with arbitrary neu- 
ron dynamics and arbitrary interactions. Employing this analysis we have investigated networks 
of IF neurons interconnected through uniform chemical synapses. In the analysis of one-cluster 
state, we have found the change of stability of a cluster, which has elucidated that in-phase syn- 
chronization of multiple IF neurons occurs only with inhibitory couplings (Fig.[2j. It must be 
noted that this analytical result well explains the structure of interneurons in the real nervous sys- 
tem, where interneurons are interconnected through inhibitory chemical synapses. In addition, 
we have investigated the entrainment of two clusters of IF neurons with different excitability 
(Fig-EJ- We have explained the narrow range of 7 e xt,2 with short decay time constant Ti in Fig. [6] 
by the bifurcation diagram described in Fig. |4] Furthermore, we have investigated one-cluster 
state of HH neurons. HH neurons show stable in-phase synchronization in the large parameter 
region of inhibitory chemical synapses, though the synchronization condition of HH neurons 
with fast decaying synapses is more complicated than IF neurons (Fig.0. 

Although van Vreeswijk et al. have proposed another type of stability criterion based on 
function G(0) [14, 15], this stability criterion is unsound in some conditions. One counterex- 
ample of their criterion is a couple of IF neurons with couplings J\\ = J22 = —J\z = —J21 = 
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g/2. With Ti = 3.5 and T2 = O.lTi, in-phase synchronization of these neurons becomes un- 
stable beyond the critical point g = 1.11 as shown in Fig. [8] While our analysis based on 
linear stability precisely yields this critical point, van Vreeswijk's criterion, namely, G(<j>) = 
-g(e- T /2) /J e Te (S[T(0 + 0)] - S[T(6 - <j))])dG with T = log (v - v r /v ), fails to give the 
critical point. Gerstner et al. have also investigated networks of IF neurons [16]. Their analysis, 
however, cannot treat the realistic form of synaptic electric current S(t) that exerts the long-time 
influence after activation within the finite size of matrix. 

The present decomposition of linear stability is simple enough to investigate the general neu- 
ron dynamics including FN neurons and HH neurons. Even when the behavior of neurons are 
chaotic [17], we are still able to evaluate the stability of cluster state using tangential Lyapunov 
exponents and transversal Lyapunov exponents [18], and such technique may give a deeper un- 
derstanding of the complicated behavior of HH neurons around the arrow in Fig. [7] It is inter- 
esting to apply the present analysis to networks including pyramidal neurons as well as interneu- 
rons [19]. The surface of the neocortex is subdivided into numerous columnar organizations, 
each of which is composed of several layers of neurons [20]. The internal and external dynamics 
of such columnar organizations would also be the future target of the present analysis. 

A The Hodgkin-Huxley (HH) equations 

The HH equations are the four-dimensional ordinary differential equations, which describe the 
spike generation of the squid's giant axon [13]. The dynamics of a neuron state vector x = 
(v, wi , W2, W3 ) T for a HH neuron is expressed as 

Cm V 
W] 

w 2 

Vi'3 

with 

«i = 
Pi = 

«2 = 

h - 
0=3 = 

ft = 

where v Na = 50 [mV], v K = -77 [mV], v L = -54.4 [mV], g Na = 120 [mS/cm 2 ], g K = 36 [mS/cm 2 ], 
0.3 [mS/cm 2 ], and C m — 1 [/xF/cm 2 ]. In the present study we set 7 ex i; = 10 ()J.A/cm 2 ) to induce 
intrinsic firing of a HH neuron. 
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Figure 1: A schematic figure explaining the definition of 8v^[t*(k)] and 8v^[t*(k)]. Membrane 
potential of a IF neuron v,- changes discontinuously at spike timing t*{k) + 8tj(k). When t *(k) + 
8ti(k) < t*(k), we define 8vf[t*(k)] by extending the solution as shown in the figure, and we 
define 8vj[t*{k)] = <5v,[f *(£:)]. When t*(k) < t*(k) + 8n{k), we define these variables in the 
opposite way. 
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Figure 2: (a) Absolute values of A.J , \\ 



and Af| for one-cluster state (Q = 1) of networks of IF 



neurons are plotted as a function of g for 7 ext j = 0,Ti = 3.5, and T2 = O.lTi. Xf always takes 

one while |Aj | is always less than one. |Afj | is less than one only when synapses are inhibitory 
(g < 0). These eigenvalues behave in the same manner even with the other decay time Ti > 0. 
(b) Z\-g phase diagram, where we fix T2 = O.lTi. A self-coupled single neuron (N = 1) has the 
stable periodic solution below g = 1. However, synchronization of multiple neurons (N > 1) 
occurs with only inhibitory couplings g < since the stability of a cluster is fulfilled with only 
inhibitory couplings g < 0. Beyond g = 1, an excessive amount of positive synaptic electric 
current leads bursting of neurons. 
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Figure 3: The result of numerical simulations with jV = 100, Ti = 3.5, and T2 = O.lTi. Dots 
represent spike timing of neurons in a stationary state, which is realized after a long run of 
simulation, (a) With inhibitory synapses g = —0.5, the perfect in-phase synchronization occurs, 
(b) With excitatory synapses g = 0.5, we observe asynchronous state, in which neurons fire 
periodically with uniformly distributed phase shifts. 
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Figure 4: Ti — q>2 bifurcation diagram for two-cluster state, where variable (fa denotes tilT 
(Q = 2, 1 < Ni = N2, g = -3, T2 = O.lTi, and / ext i = / extj2 = 0.) The solutions represented by 
thick lines satisfy the stability of mean state and stabilities of clusters, while solutions represented 
by the dotted lines lack one or both of stabilities. The out-of-phase solutions plotted by the 
thin dotted line (Tj < 2.8) is invalid since in these solutions membrane potential v, crosses the 
threshold 9 multiple times. 




Figure 5: Entrainment of two clusters of neurons with different excitability. The solution <j!>2 is 
plotted against / extj 2 for the fixed value of 7 ext ,i = (Q = 2, 1 < Ni — N2,g — — 3, Ti =3.5, and 
T 2 = 0.lTi.) 
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Figure 7: T\ —g phase diagram for one-cluster state (Q — 1) of multiple Hodgkin-Huxley (HH) 
neurons (2 < AO under the condition T2 = 0. 1 Ti . "s" ("u") in the figure indicates the region for 
the stable (unstable) one-cluster state. Around the arrow we find a lot of isolated regions for the 
stable one-cluster state. Note that we apply constant external electric current / ext = 10 (jxA /cm 2 ) 
to all of HH neurons so as to induce intrinsic firing of neurons. 
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Figure 8: Stability of in-phase synchronization of a couple of IF neurons interconnected with 
■^n = J22 = —J\2 = —J21 = g/2 (Ti = 3.5 and T2 = O.lTi). (a) Absolute values of kf and Aj are 
plotted as a function of g. Beyond g = 1 . 1 1, the in-phase synchronization becomes unstable, (b) 
The result of numerical simulations with g = 1 .0. A couple of neurons show in-phase synchro- 
nization, (c) The result of numerical simulations with g = 1 .2. Only a single neuron fires at high 
frequency in the winner-take- all fashion. 
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